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. A non-perturbative algebraic theory of lattice Boltzmann method is developed based on a symme- 

try of a product. It involves three steps: (i) Derivation of admissible lattices in one spatial dimension 
through a matching condition which imposes restricted extension of higher-order Gaussian moments, 
CNj , (ii) Special quasi-equilibrium distribution function found analytically in closed form on the product- 

lattice in two and three spatial dimensions, and which proves factorization of quasi-equilibrium 
moments, and (iii) Algebraic method of pruning based on a one-into-one relation between groups 
of discrete velocities and moments. Two routes of constructing lattice Boltzmann equilibria are 

■ distinguished. Present theory includes previously known limiting and special cases of lattices, and 
f"^ ' enables automated derivation of lattice Boltzmann models from two-dimensional tables, by finding 
£f} , roots of one polynomial and solving a few linear systems. 

PACS numbers: 47.ll.-j, 05.20.Dd 
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I. INTRODUCTION 

g; 

There were a few recent attempts [l|, 0J!, [3, 0, [f| @] to construct a theory of the lattice Boltzmann (LB) method 
- a modern approach to fluid dynamics [8|. This is due, in the first place, because LB models currently in use are 
not "sufficiently" Galilean invariant (the feature that LB improved on from its predecessor, the lattice gas model, but 
failed to resolve completely). Even though the Galilean non-invariance of current LB models was very well known 
L | right from the beginning @, curing this drawback resisted for a long time. Insufficient Galilean invariance of the 
very basic LB at a constant temperature is a precursor of many difficulties, in particular, i n ap plications of LB to 
Q high Reynolds number hydrodynamics [HI, [12] , multi-phase flows [l| and compressible flows Q EH] . It is quite well 
O , understood that the current "standard" LB models are too much constrained by the "small" number of the discrete 
velocities, and lattices with "more" velocities are required in order to overcome these limitations. However, early 
I \ attempts to introduce lattices with more velocities were unsuccessful because of a severe numerical instabilities of the 

■ resulting LB schemes [MElEEI- 

^\ ' Important progress was recently achieved in [J Q , where the construction of the higher-order LB was formulated 
. as the construction of the entropy In particular, [H, Q explained why some of the most obvious suggestions for 
' higher-order lattices are bound to failure due to the fact that no entropy can be constructed for them. The entropy 
construction of Refs. p], Q has led to admissible lattices in three dimensions which enable LB models with better 
t-H ■ properties but derivation of such lattices (the procedure termed pruning in Ref. Q ) remained a rather tedious search 
among large families of lattices. Apparently, some kind of simplicity was still missing at that stage, and a fully 
analytic approach to pruning is a challenging task. On the other hand, symmetry with respect to a group of rotations 
was invoked recently for a classification of isotropy of higher-order LB models (in two dimensions) [201 ] - However, 
the information about the isotropy of the higher-order lattices alone is insufficient if we want to address stability (or 
instability) and the form of the equilibrium on each specific lattice. 

In this paper, we develop a theory of higher-order LB based on a symmetry of a product. We remind that such a 
symmetry is deeply rooted in the classical kinetic theory since its beginning, the seminal Maxwell's derivation of the 
equilibrium of the three-dimensional. Isotropy (independence of the equilibrium on the direction) in Maxwell's famous 
derivation comes from the fact that the product of one-dimensional Maxwell distributions depends only on the isotropic 
quantity, the kinetic energy of the particles: exp(— v%) exp(— Vy) exp(— v\) = exp(— v ■ v). Our consideration of the 
lattice Boltzmann method is based on theproducts of one-dimensional functions. The present theory of LB method 
is algebraic (rather than group-theoretic [20. [2l|or function-theoretic P, 0) and non-perturbative (it is not based on 
polynomial expansions of the Maxwellian 0, 0j S S 0] ) • The latter is important for preserving the symmetry of the 
product, as we will see it below. The resulting theory is remarkably constructive and simple, and consists of three 
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major steps: The construction begins in one dimension where we identify admissible lattices (sec.|TI|. At this step, we 
reveal the reference temperature (of the Maxwell distribution represented by the given one-dimensional lattice). This 
information is then immediately transferred (sec. Illlj) into three dimensions with the help of a special unidirectional 
quasi-equilibrium on a "large" lattice formed by all possible direct products of one-dimensional velocities (Maxwell 
lattice) (for general issues related to quasi-equilibria see [22| )■ The result of sec. IIIII (see Eq. Q below) extends the 
product form onto the entire quasi-equilibrium populations. The advantage of the unidirectional quasi-equilibrium 
on product-lattices is twofold: It has a simple structure of the corresponding moment representation (see Eq. ([5]) 
below) , and constructing the equilibrium is a mere substitution of the one-dimensional data for one-dimensional non- 
conserved moments. We distinguish between two routes to obtain the equilibrium for lattice Boltzmann models: The 
equilibration (minimization of the entropy function under constraints of local conservation) and the Maxwellization 
(promotion of Maxwell's equilibrium values for the non-conserved moments). The Maxwell lattice is an "ideal" lattice 
in three dimensions, it replicates all the information gained in one dimension. "Ideal" also means that the information 
about three dimensions is represented without correlations in the product-form ((4]). Based on the results of sec. IIIII 
in sec. |IV] the analytical method of pruning is developed. The main ingredient in this approach to pruning is the 
two-dimensional key-table which furnishes the one-into-one relation between groups of velocities and moments, and 
which is relatively easy to analyze even for large velocity sets. The pruning algorithm is explained with the examples 
of the familiar D3Q27 lattice and the higher-order D3Q125 lattice. In particular, the Maxwellization based on the 
pruning of the unidirectional quasi-equilibrium moment system, derives equilibrium distributions by solving linear 
algebraic systems. Finally, the results are discussed in sec. |V] 



II. MAXWELL LATTICES IN ONE DIMENSION 



Since our construction will be based on the one-dimensional lattices, it is important to sort it out right from the 
beginning which one-dimensional velocity sets are admissible, and which have to be rejected. Therefore, we consider 
the one-dimensional sets of discrete velocities V, with Q the total number of the velocities (below, we consider Q odd 
but same considerations apply also to Q even) . The discrete velocities Vu\ € V are assumed integer- valued such that 
V(i) = i. The basic mirror symmetry of V assumes that if Vu\ € V then also i>(-j) <E V, and thus stopped particles 
with t>(o) — arG always included. Corresponding populations are denoted f(j\, and we use convenient normalization, 
f(i) = PVU)- Summation over discrete or integration over continuous velocities will be denoted as (. . . ), thus p = (f(i))- 

Discrete velocities V are so chosen as to reproduce the moments of the one-dimensional Maxwell distribution 
function, 



where 
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1 (normalization), 
u (flow velocity), 

To + u 2 = II M (equilibrium pressure at unit density), 



3T u + u 3 = 
3T 2 + 6T u 2 



q (equilibrium energy flux at unit density), 



R 



M 



l5T 2 u 



WTqu 3 



and so on, to which we refer as Maxwell's (M) moment relations. 

The first information revealed from the lattice is the reference temperature To at which (a part of the) Maxwell's 
moment relations will be verified. This is done with the help of the closure relation and the matching condition. 
The closure relation for the set V with Q velocities (Q odd) is a linear relation between the Q-th power of the 
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Closure 
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3 


{0, ±1} 


(!) W 


1/3 


5 


{0,±1,±3} 


v fi) = lOwffl - 9«(i) 


1 ± V 2 /5 


7 


{0,±1,±2,±3} 


ww = Uv h - 4 H> + 36 «w 


0.697953 


9 


{0,±1,±2,±3,±5} 


wfo = 39ufo - 399ufo + 1261ufo - 900« w 


0.756081, 2.175382 


11 


{0,±1,±2,±3,±4±5} 


= 55w^ - 1023u f ' o + 7645vfo - 21076vf o + 14400u w 


1.062794 



TABLE I: One-dimensional Maxwell lattices with odd number of integer-valued velocities, Q = 3,5,7,9,11. Second column: 
Lattice vectors; Third column: Closure relation, defining the reference temperature To through the matching condition (fourth 
column) . 

velocities, uj^ , and the lower-order odd powers, starting with and ending with . Such a linear relation 

always exists, and reflects the fact that only Q velocity polynomials, . . . tW^ 1 arc linearly independent. For 

example, for V = {0, ±1} (D1Q3), the closure relation is vf = Urn (cube of any velocity from the D1Q3 set is the 
velocity itself), for V = {0, ±1 ± 3} (D1Q5) it is vf-\ — 10u?^ — 9v^, and so on. The existence of the closure relation 
implies that the moment M(q) = (ip^yv^) cannot be assigned at one's will, and that only the linear in u term of 
this moment at equilibrium can be made consistent with the corresponding Maxwell's value M^. This leads to 
the matching condition which decides about the reference temperature T . For example, for D1Q3 the third-order 
moment M( 3 ) = {ni) v fi)) equals M( 3 ) = u for any population set, equilibrium or not. On the other hand, the 
Maxwell's expression, M^£ = 3T u + it 3 , contains also the cubic term it 3 which cannot be made consistent with 
the previous expression. Only the linear term can be matched, 3T u = u, if the reference temperature is set to 
To = 1/3. Similarly, for D1Q5, = \bpT^u + 0(u 3 ), thus, the matching condition for linear terms becomes 

15TQ — 30T + 9 = 0. The latter equation reveals two values of the reference temperature, T = 1 ± */2/5. This 
example also explains why the shortest D1Q5 lattice is {0, ±1,±3} and not {0,±1,±2}: For the latter, the closure 
relation is ujL = 5u?^ — 4i>( i ) , and T is found as a solution of 15Tq — 15T + 4 = which has no real- valued roots, and 
hence does not define any reference temperature. This procedure is immediately applicable to any lattice (Appendix 
1a"]) . In Tablc[H wc collected Maxwell lattices with Q = 3,5, 7, 9, 11, together with the corresponding closure relations 
and the reference temperatures. 

Once the reference temperature is revealed, we immediately derive the equilibrium values of the populations at 
u = and unit density (weights) Wuy 

For this, we introduce the complete set of moments (at unit density): 

M (0 ) = 1 = M(i) = (<P(i)t>(i)), ■ ■ ■ , M( Q _ 1) = (<P(i)i$) -1) ). 

Denote M. — {Mm, . . . ,Miq_\\\ the totality of the moments, excluding M(o) = 1, and f2 the set of their values at 
which the solution to the latter Q x Q linear system is positive. This solution <P(i)(M) is always easily found from 
the above Q x Q linear system, and we denote 

f(i) =PV(i)(M). 

For example, for D1Q3, functions iffj\(u,Tl) are found as the solution to a 3 x 3 linear system, tpr ^ — (1 — II), 
ip {±1) = (1/2) (Hi it), where u = M (1) and IT = M (2 ), while O = {u, II : < II < 1, |u| < IT}. For the D1Q5, we need 
two more moments of order three and four, q = (f^v^y) and R = (ip^yv^y), and thus 

no) = i-n + i(i?-n), 

n±D = ^[M9u-q)-R + 9H\, 

±(q-u) + ^(R-IL) , 

and so forth. In order to reveal the weights, we substitute the equilibrium values of the corresponding moments at 
u = into the above formulas for a ^ -^o already available, to derive 
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This gives W(p) = 2/3, W {±1) = 1/6 for D1Q3, W {0) = (4/45) (4 + y/T5) , W (±l) = (3/80) (8-\/I6), W c±3 ) = 
(1/720) (16 - 5-y/TO) for D1Q5 (at T = 1 - ^2/5) and so on. 

Once the weights and the reference temperature are derived, we can immediately proceed with the evaluation of 
the equilibrium populations. There are two options: 

(i) Equilibration. The weights Wu) > define the entropy function H = (fu,) In^f^/Wu))}. The equilibrium 
populations /Sx = pf E ^ are defined as the minimum of H, conditioned by density p and velocity u. Let us 
distinguish between the velocity u and the higher-order moments by writing 

M = {u,Af}, 

so that 

Af = {M [2) ,...,M (Q _ 1) }. 
The above functions tpu\(u,J\f) are substituted into H to give H(p, u,Af) = p\np + pH(u,Af), where 

H(u,N) = <v W («,A0ln (¥>(i)(u,JVW«))> • 
The equilibrium is found from the equations, 

dH{u,N) dH(u,Af) 
8M m ~ ' '"' <9M(q_x) ~ 

These equations define the equilibrium solution A/" E = 7V E (u) . Exact solution is available (so far) only for the 
D1Q3, where Af consists of the pressure II only; then IT 3 = (l/3)(2yT+~3l4 _ i). i n ot her cases, various 
solution procedures can be readily applied to get approximations to M 0] . The equilibrium is thus 

f^=pip (l) (u,Af E ), 

where AT E is exact or approximate solution to the extremum condition. 

(ii) Maxwellization. Alternatively, we can promote Maxwell's expressions of the moments AT m (Tq,u) to derive a 
different set of equilibrium populations, 

For example, the Maxwellization of the D1Q5 model is accomplished upon substitution of II M , q M and R M in 
the above expressions for <f(o), V(±i) an d <P(±3) : 

/(^ = P {l - (To + u 2 ) + l [(3T 2 + 6T u 2 + u*) (To + u 2 )} | , 
/(±i) = J^P {±[9« - ( 3T o« + u 3 )} - (3T 2 + 6T a u 2 + u 4 ) + 9(T + u 2 )} , 
/(±3) = |±[(3T oM + u 3 ) -u} + i[(3T 2 + 6T u 2 + u 4 ) - (T + M 2 )]| , 
with T = 1 - ^275. 

Maxwellization is easier than equilibration since functions Af M are known from the one-dimensional Maxwellian f^f. 
In a contrast to the continuous velocity case where f E = f™, Maxwellization is not the same as the equilibration on 
the lattice. In order to illustrate this point, we present a comparison of the lattice Bhatnagar-Gross-Krook (LBGK) 
simulation of a one-dimensional shock propagation with two different equilibria. In Rcf. [lj, it was shown that the 
LBGK model on the Maxwell D1Q5 lattice V = {0,±1,±3}, with the equilibrium constructed by the equilibration 
procedure (that is, via the entropy minimization ) is superior in terms of numerical stability to the LBGK model on 
the inadmissible lattice V = {0, ±1, ±2} of Ref. [Tj|. In Fig. [TJ we present the result of the same simulation for the 
LBGK model with the equilibrium obtained by Maxwellization. Both simulations, with / E [l| and / M (present) agree 
well with each other, and show the same stability properties. 

Thus, the one-dimensional decoding is complete, we have derived reference temperatures and weights for an arbitrary 
onc-dimcnsional velocity set just from the lattice itself. In the next step we are going to transmit the one-dimensional 
information into three dimensions. We close this section with a few comments: 
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FIG. 1: (Color online) Shock tube test: Comparison of the Maxwellization on the D1Q5 lattice V = {0, ±1, ±3}) (line) with 
the equilibration of Ref. [IJ on the same lattice (symbol). Initial condition for the simulation was a density step, p = 3.0 
for x < L/2 (L being the length of domain), p = 1.0 for x > L/2 (same as in [l], 01|). The snapshot of the density profile 
corresponds to kinematic viscosity u — 0.138 (the LBGK model of Ref. [l8| on the inadmissible lattice V — {0, ±1,±2} is 
unstable at this viscosity and is not shown). 

• The reference temperatures for the Maxwell lattices collected in Table [H and the corresponding weights, coincide 
with those found in Refs. [2, [3] using the entropy construction [l9j |. The entropy construction derives the weights 
and the reference temperatures by comparing higher powers of velocity of A4 E to Maxwell moments. However, 
the present derivation via closure relation and matching condition is more direct and simpler. While the 
coincidence of the results obtained by two methods is quite remarkable, and suggests that the two approaches 
may be equivalent, the full proof of this statement is not available at the time of this writing, and is left for a 
further study. 

• It should be stressed that the three- velocity case (the basis of the "standard" LB models) is an exception: any 
set V = {0, ±r} is Maxwellian (the corresponding closure relation, v%*. = r 2 v^, results - through the matching 

condition - in only a trivial re-scaling of the reference temperature, Tq = r 2 /3). With Q > 3, by far not 
every lattice is Maxwellian (see the above example of V — {0, ±1,±2} and Appendix IB"]) . The above concept 
of Maxwellization applies exclusively to the Maxwell lattices. Additional comments on Maxwell lattices and 
matching condition will be given in sec. |V] 

• The two values for the reference temperature for Q = 5, 9 (Table [TJ) correspond to a Gaussian-like shape of the 
weights (Wu\ < Wyi if \i\ > \j\) for smaller To, and to a non-Gaussian shape for larger To (cf. Ref. Q ). Below, 
we consider To corresponding to the Gaussian-like case in all the examples. 



In three dimensions, we first construct the product-lattice (or Maxwell lattice), induced by the one-dimensional 
Maxwell velocity set V, that is, 

(i) The velocities are direct products of one-dimensional velocities, 



III. MAXWELL LATTICES IN THREE DIMENSIONS AND UNIDIRECTIONAL 

QUASI-EQUILIBRIUM 



Unidirectional Quasi-Equilibrium 



V(i,j, k ) = («(i).«Ci). «(*)), 
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(ii) Corresponding weights are algebraic products of the one-dimensional weights, 

WW) = W (i) W {j) W {k) . 

The entropy on the product-lattices is defined as 

// (j,,, h.f„. f 4r w ))■ 

\ \W(i)W (j) W (k) J / 

Moreover, the density is defined in the usual way, p = (f(ij.k)): ancl we introduce a set of special unidirectional 
moments M a defined as 

pM x[n) = {f(i,j,k)V^)), pM y[n) = (/(i,j,fe)«y)}, pM z[n) = n = 1, . . . , Q - 1. (2) 

Note that M a ^ = u a are the components of the three-dimensional velocity, while the rest of the unidirectional 
moments are the diagonal components of the corresponding tensors. For instance, M Q ( 2 ) are the diagonal components 
of the pressure tensor (at unit density), M a ^ - of the third-order moment tensor etc. 

Using ([2]), we define a special unidirectional quasi-equilibrium state (UniQuE) as the minimizer of the entropy 
function under the constraints imposed by fixed density and fixed unidirectional moments That is, UniQuE 
populations f*(p,M x ,M y ,A4 z ) are defined as the solution to the variational problem, 

H -> min, (f^k)) = P, = P M x(n), (/(i,j,fc) u ^)) = P M y(n)> (f(i,j,k)v"k)) = pM z(n) , n = 1, . . . ,Q- 1. (3) 

The central result of this section is given by the following Theorem: 

Solution to the conditional minimization problem (|3|) is explicitly given by the formula 

f(i,j,k) = P<P(i)(M x )(p u) (M y )(p {k) (M z ), (4) 

where the positive one-dimensional populations Lp^{M a ) arc defined by solving the one-dimensional moment system. 

To prove this (see Appendix[C]), it is sufficient to notice that the solution to the minimization problem in terms of the 
Lagrange multipliers reduces to three decoupled one-dimensional problems of the form, (<^(i)) = 1, {^(1)^)) = -^x(n)) 
and similarly for y, z. Solution of each of these problems is given by the unidirectional functions tp(A4 a ) discussed in 
sec. El 

UniQuE (J3J) is a family of populations defined by 3Q — 2 parameters in the Q 3 -dimensional space, and is a fully 
factorized population: In order to construct (|3]), we plug M. a instead of M. in the one-dimensional functions (p^)(A4), 
and multiply results for various i and a. 

The above theorem about UniQuE applies to any Maxwell lattice. We note in passing that special versions of 
UniQuE for the two-dimensional D2Q9 lattice was constructed in [l5[ and [2^|, and for the D3Q27 - in [24] from 
the direct minimization of entropy. UniQuE ^ is the most crucial element in passing the information to three 
dimensions. Note that, in general, it is impossible to find closed-form expressions for a quasi-equilibrium which 
minimizes the entropy under arbitrary constraints. UniQuE is the exceptional case because the solution is induced 
by the one-dimensional solutions which are explicitly known. This is possible only with the special choice of the 
constraints (unidirectional moments), and only on Maxwell lattices. This has a few immediate implications, two of 
which will be mentioned now. 



Moment representation 

The product-lattice generated by Q one-dimensional velocity vectors is characterized by Q linearly independent 
moments, 

pMi mn = (/(jj.fc)^)^™)^)), l,m,n G {0,...,Q- 1}. 

On the other hand, UniQuE is fully described by only 3Q — 2 moments (density and unidirectional moments M;oo = 
M x ([) etc). Thus, the rest of the moments become functions of density and unidirectional moments when evaluated 
on the UniQuE Evaluation is straightforward thanks to the product-form of the latter: 

M hnn = M x (l)M y ( m )M z ( m ). (5) 

Thus, the moment representation of UniQuE ([5]) is a simple algebraic rule: One considers all possible products of 
functions M a ( p ) with different spatial index a, times the density p, where the number of functions in each such product 
does not exceed three. Example of the UniQuE moment system ([5]) for D3Q27 Maxwell lattice is given below in Table 
IIIII Finally, since the moment and the population representations are equivalent to each other, we can now read ([5]) 
"from the right to the left" and say that it defines UniQuE upon inverting the Q 3 x Q 3 linear system ((5]) with the 
specified right hand side. This remark will be important later when we will consider sub-lattices of the product-lattice. 
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Equilibration 

The term " quasi-equilibrium" in the notion of UniQuE means that it is "less equilibrated" than the equilibrium. 
The equilibrium (at the fixed reference temperature) minimizes entropy under fixed density and velocity u a = M a (jy 
Let us distinguish between the velocity u a and the higher-order moments by writing M. a = {u a ,J\f a }. The above 
theorem about UniQuE implies that the following two routes to equilibrium are equivalent: 

• The direct equilibration through minimization of H under fixed p and u a , and 

• The two-step equilibration, of which the first step is the " quasi-cquilibration" by minimizing H under fixed p 
and M. a (resulting in UniQuE), followed by the second equilibration step during which the UniQuE entropy 
H* = plnp + p[H(u x ,Af x ) + H(u y ,J\fy) + H(u z ,J\f z )] is minimized with respect to Af a under fixed p and u a , 
a = x,y, z. 

It is obvious from the product-form of UniQuE ^ that the second minimization reduces to the one-dimensional 
equilibration of sec. [12 and thus 

/(L'.fc) = W(0(«*.^)*'M(«v.^^)K,^). (6) 

This again requires only the input from the one-dimensional lattice (functions Af^j)- In other words, the UniQuE 
becomes equilibrium when the equilibrium values of the unidirectional moments are substituted into (|4]). A few 
comments are in order: The lattice Boltzmann equilibria on the product-lattices are constructed in such a way that 
the higher-order tensorial moments of a certain order render isotropic (to a certain order in the powers of the velocity 
components u a ) Q. On the contrary, the special quasi-equilibria considered above are anisotropic (their construction 
is based explicitly on a fixed Cartesian system of coordinates which is manifest in our choice of the parameters, 
the unidirectional moments). Yet, the evaluation of these anisotropic moments (which are typically the diagonal 
components of the corresponding higher-order tensors) at the equilibrium renders the same degree of isotropy for the 
entire moment tensors at the equilibrium. Or, in other words, the control (bringing to the equilibrium) over just 
the diagonal components of moment tensors is sufficient to control the entire tensors (including various off-diagonal 
components which are not explicitly targeted in the construction of the quasi-equilibrium) . This fully corresponds to 
Maxwell's argument on how the equilibrium in the three-dimensional gas become isotropic based on the independence 
of the three directions. 



Maxwellization 

Same as in sec. [TTl there is a different route to define the equilibrium on the product lattice by simply plugging in 
Maxwell's values AT M (To,u) into UniQuE to get a three-dimensional Maxwellization, 

/&.*) = P^(u x ,^) m (u y ,^) m (u z ,^). (7) 

Note that ([7]) is not the same as (H2). Moreover, ([7]) differs also from the standard polynomial equilibrium on the 
product-lattices (for example, for the D3Q27, is a polynomial of the order six, while it is a second-order polynomial 
in the standard LB model). As an illustration, we collected all the populations mentioned so far (UniQuE, equilibration 
and Maxwellization) for the D3Q27 in Appendix iDl 

Discussion 

Thus, the transmission of the one-dimensional information to three dimensions is now completed for the Maxwell 
lattice. Arguably, this is a transmission "without errors", all the information about the Maxwell's relations collected 
for the one-dimensional distribution is manifest in the three dimensions once the product-lattice is used. For example, 
the Maxwellization on the Maxwell lattices (J7J recovers Q 3 moments of the three-dimensional Maxwellian: 

MZn = Mf {l) M^ m) Mf {m) . (8) 

Moment relations (J5J) set the maximal possible accuracy achievable on the Maxwell lattice (for example, the moment 
system as recovered by the kinetic equation dtf + v ■ V/ = — (l/r)(/ — / M ) is the closest approximate to a truncated 
moment equations system of the Boltzmann equation with the Bhatnagar-Gross-Krook collision operator). Note 
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TABLE II: Backbone moments of the D3Q125 UniQuE system (f5J arranged in columns according to their order. Upper left 
corner displays the backbone moments of the D3Q27 (see Eq. ([9])). 

that, in general, product-form of equilibria such as ((7J) or (J6j> should be preferred in LB computations [25|. Whereas 
LB equilibria found by other methods (in particular, those using a polynomial expansion of the Maxwellian and 
quadrature approximations [H Hi HI Q ) can be recovered upon a further expansion and neglect of higher-order terms 
in ([7]) or ([S]), this discussion remains out of scope of the present paper since the method used here unambiguously 
results in the product-forms ([7]) and 

The drawback, however, is that the number of the velocities needed for this "error- free" transmission grows as Q 3 
which becomes a large number. Therefore, we need to consider an "incomplete" transmission by sacrificing some of 
the moments and reducing the number of velocities accordingly (pruning). Above, we have remarked that UniQuE of 
the product-lattice can be computed from the full Q 3 x Q 3 linear moment relations (f5|). However, if we consider a part 
of the moment system ([5]) including moments of primary importance to the hydrodynamics only, this computation 
can be accomplished with a lesser number of the populations, or, equivalently, with a lesser number of the lattice 
velocities. In view of a large number of different moments, how to do this in a systematic fashion? The answer to 
this question is central to the present theory, and will be given in the next section. 

IV. PRUNING AND SUB-MAXWELL LATTICES 
A. Backbone moments and sub-Maxwell lattices from Key- Table 

Looking back at ([5]), we notice that only even-order moments give a non- vanishing contribution to this system at 
the equilibrium at velocity equal to zero. Indeed, the odd-order moments such as u ai Q„a„ or Q^ 7 , etc. all vanish 
at u a = 0. What remains are the even-order moments which we call the backbone moments. These are various 
even-order unidirectional moments and various products constructed with their help, up to the triple product of the 
highest-order even unidirectional moments. For example, for the D3Q125 lattice (the Maxwell lattice generated by 
the one-dimensional velocity set V — {0, ±1, ±3}) there are ten different types of the backbone moments arranged in 
the increasing order from zero to twelve (see Table HI)) . On the other hand, the product lattice can be represented 
as a collection of shells, each shell contains all the velocities with the same magnitude and symmetry with respect to 
reflections at the origin and permutation of components. It is important to realize that 

The number of different types of the backbone moments equals the number of shells. 

This observation makes it possible to find a one-into-one relation between the backbone moments and the velocity 
shells which has a form of a two-dimensional key-table (KT) . Let us explain its construction with the example of the 
D3Q27 product-lattice (see Eq. ©). 

The backbone moments are then of the four types: 

i, n Q , n a n^ (a^/3), n x n y n z . 

At the zero- velocity equilibrium / oq (where eq is either E or M), these are four different values, 

i, ic = To, rcrr^ = t 2 (a ± (3), n^n^n^ = t 3 . 

On the other hand, the D3Q27 lattice is composed of four shells: 

V Q = {(0,0,0)}, 

V 1 = {(±1,0,0), (0,±1,0),(0,0,±1)}, 

V 2 = {(±1,±1,0),(±1,0,±1),(0,±1,±1)}, 

V 3 = {(±1,±1,±1)}. 

The shells V s enumerate the rows in the KT ©. Now, we compute contribution of each shell to each backbone 
moment, introducing the (yet) unknown weights W s for the velocities of each shell. This corresponds to the 4x4 
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entries of KT (|9|). Next, summing up the entries in each of the four columns, and equating the result to the equilibrium 
value of the corresponding moment, we get a 4 x 4 linear system for the weights, Wq + QW\ + 12W 2 + 8W3 = 1, 
2Wi + 8W 2 + 8W3 = ?o, AW 2 + 8W 3 = T 2 , 8W 3 = T 3 . This system is what remains from © of the D3Q27 at 
zero-velocity equilibrium. Substituting T = 1/3, we get W = 8/27, W\ = 2/27, W 2 = 1/54, W 3 = 1/216, the result 
which we already knew from the product-form. 



•S 


1 


n a 


n a Uf3 


(a 







W 













1 


6Wx 


2Wi 










2 


12W 2 


8W 2 


4W 2 







•3 


8W 3 


8W 3 


8W 3 




8W 3 



(9) 



Now, the pruning method with the help of KT ([9]) consists of erasing one or several rows and of the same number of 
columns. Erasing rows is the pruning of the lattice by shell wise discarding of the velocities, whereas erasing rows is 
sacrificing some of the backbone moments, that is, reducing the accuracy of the UniQuE moment system (|5|). Lattices 
constructed in this way from a Maxwell lattice will be termed sub-Maxwell lattices. Certainly, in order this procedure 
to be useful for a further construction of LB models, we should favor lower-order moments as they contain most of the 
information about the hydrodynamics. In the present illustrative example of the D3Q27 product-lattice, this means 
that we should keep the first and the second columns (corresponding to the density and to the diagonal components 
of the pressure tensor) since these are required for recovering the Navier-Stokes equations at low Mach numbers, while 
the higher-order moments (last two columns) can be sacrificed in the pruning procedure. 

It is easy to see how the "standard" LB lattices come out as the result of this process. Erasing the last row (s = 3) 
and the last column (rLriyllz) in KT ([9]), summing up the remaining columns, equating the results to the values 
of the corresponding backbone moments at zero velocity equilibrium at the reference temperature To = 1/3, and 
solving the resulting 3x3 linear system, gives the weights Wo = 1/3, W\ = 1/18 and W 2 = 1/36 which describe 
the "standard" D3Q19 lattice. Erasing the third row (s = 2) and again the last column gives Wo = 2/9, W\ = 1/9 
and W 3 = 1/72, which is another standard D3Q15 lattice. Finally, a less standard D3Q13 lattice [2y] corresponds to 
erasing the second and the last rows (s = 1 and s = 3), and two columns, next to the last and the last (II Q n^ and 
n^nyllz), resulting in Wo = 1/2, W2 = 1/24. Note that the present examples illustrates a complete pruning: the 
three lattices just mentioned are the only pertinent to recovering the isothermal Navier-Stokes equations as the result 
of pruning of the D2Q27. 

Key-tables similar to ([9]) are obtained in a straightforward manner for product-lattices with any Q, and are relatively 
easy to analyze (for example, for the D3Q125 lattice, the number of types of the backbone moments is an order of 
magnitude less than the total number of moments, cf. Table ITT| . Following this procedure, we easily identify, for 
example, the recently introduced D3Q41 lattice @: The six types of the backbone moments retained are: 1, II Q , 
n Q rLi, R a , H a R/3 and n T II y rL. The retained six shells include the four shells Vq, . . . , V 3 of the D3Q27 mentioned 
above together with V4 = {(±3, 0, 0), (0, ±3, 0), (0, 0, ±3)} and V 5 = {(±3, ±3, ±3)}. The two latter shells contain 14 
velocities which, added to the 27 make up the D3Q41 lattice. Computing the contribution of these six shells to the six 
backbone moments, and solving the resulting 6x6 linear system, we immediately obtain the corresponding weights, 



Wo 
Wi 

w 2 
w 3 
w± 
w 5 



1 - T [270 - T (263 + 1021b)], 

iTo[9-T (12 + 13T )], 
16 

2 °' 

^T 2 (9-19T ), 

To(3T -l)(9-T ), 



(10) 



1296 
1 

5184 



T 3 (3T - 1), 



which are positive at To = 1 — \/2/5 (see Table HJ, and coincide with those reported in 0]. 

The pruning of the Maxwell lattice using its KT derives the important information, the weights W s corresponding 
to the retained shells (for the pruned lattices, the weights are not products of any one-dimensional weights any 
longer, as it was for the product-lattice). This immediately triggers the option of equilibration by minimizing the 
corresponding entropy 0. The equilibration is performed under fixed density and velocity, which are now defined on 
the sub-Maxwell lattice. 
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TABLE III: Projection pruning of the D3Q27 UniQuE moment system (0. Moments are grouped in columns, according to 
their order, from to 6. Each row contains the moments retained by a particular lattice. Backbone moments are indicated 
first. First row (D3Q27) represents the full UniQuE moment system Filling out the rows corresponding to D3Q19, D3Q15 
and D3Q13 is explained in the text. Maxwellization (construction of the equilibrium) is achieved by replacing il Q — > E[Jf, where 
Il„ = To + u„, and To = 1/3 is the reference temperature. Example of D3Q19 is presented in Appendix iDl Eqs. (|D7[) and 
(|D8I) . Example of D3Q13 is further discussed in Appendix [El 

Finally, wc remark that KT establishes the most "fine-grained" (one-into-one) correspondence between (groups of ) 
velocities and moments (it is not possible to establish a "finer" correspondence between the moments and the velocities 
than that provided by KT since many velocities contribute to each particular moment). The relation between velocity 
shells and backbone moments, as presented by KT, is therefore the optimal setting for pruning, in general. 

B. Projection pruning 

The advantage of the above entropy pruning (EP) is that, once the weights are found from KT, we do not need to 
care about the higher-order moments since their equilibrium values will be decided by the corresponding equilibrium 
/ E . The disadvantage is that we (still) need to solve a nonlinear minimization problem to find f E . Therefore, a 
different way of pruning can be offered which avoids the entropy minimization and is much easier to execute. 

This route is, in fact, a continuation of the KT to include a part of the moment system ([5]), addressing also the 
moments which were washed out at the zero-velocity equilibrium. Let us again explain it with the example of D3Q27 
(see Table ITTTI) . First, we group all the moments §5§ according to their (usual) order from to 6 (the highest-order 
moment corresponds to the triple product n :E ITj,IT 2 ), writing the backbone moments first (first row of Table ITTT)) . 
For each lattice found from the above analysis of KT, we fill out the corresponding row by retaining (a part of) the 
moments ([5]), moving from the left to the right (from the lower to higher order moments). For example, for the 
D3Q19 lattice (second row in Table HIT)) . we first include all the moments in the columns 0, 1 and 2 as they define 
the basic fields (density and velocity), and the pressure tensor. In the column 3, we can include all the third-order 
moments except for u x u v u z because Mm degenerates on the shells retained in the D3Q19: Since any velocity vector 
of D3Q19 contains at least one zero component, we have v^v^v^ — (i 7^ j 7^ k) for any vector. This degeneracy 
precludes the moment = u x u y u z to be retained by the moment system of D3Q19, and we proceed to the next 

column, where we can retain only the three backbone moments. In the case of D3Q15, the situation is opposite at 
the column 3: while the moment Mm is non-degenerate, and thus the value M-j^ can be now retained, the three 
pairs of moments, M120 and M102, M210 and M012, and M021 and M201 become degenerated, and only the three 
linearly independent combinations can be retained. For that, we choose symmetric combinations, as shown in Table 
IIIII Finally, the three backbone moments are degenerated by D3Q15, Af 2 2o = M202 = M 02 2, and we are able to retain 
their symmetric combination. Similar considerations apply also for the last (D3Q13) lattice reported in Table ITLTI 

Now, the number of retained moments in each row of Table IIIII equals to the number of the populations of the 
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corresponding lattice. Consequently, these moment relations, with the right hand side given by Table IIII1 can be 
readily inverted to derive an analog of the UniQuE, 

f*i,j,k) = PV^ 3 ,k)(M x ,M y ,M z ), (11) 

where now k) spans not the whole range of indices but only those corresponding to the retained shells. Conse- 
quently, <pu,j k) do not have the form of a product of the unidirectional functions ((H) (although it resembles the latter, 
as illustrated by the D3Q19, see Eq. (|E)7|) in Appendix |D|) . Function /* (fTTj) represents the UniQuE moment system 
([5]) in the best possible way allowed by the reduced number of velocities, thereby providing a projection of the D3Q27 
lattice onto the corresponding pruned lattice. For that reason, we term the present method as projection pruning 
(PP), in order to distinguish it from the entropy pruning. 

Since PP derives /* (fTTj) from the moment system of the Maxwell lattice (O, the notion of the equilibrium for it is 
also a derivative of the corresponding results for the UniQuE (|T]) : It is either equilibration, induced by the equilibrium 
values Ai® — {u a ,N^} of the corresponding one-dimensional Maxwell lattice, 

/(i'.fc) = PV{i,i,k) (K,^}, {UyMy}, {UzMf}) : (12) 

or Maxwellization, induced by the Maxwell values of the same one-dimensional moments A4„ = {)i a ,A^} 

fc=W(W) (K,AA x M },K,AA y M },K,A^}). (13) 

In Appendix HO we give example of /* (HU and f M (UJ) for the D3Q19 sub-Maxwell lattice (Eqs. (jDZj) and (|D8|) . 
respectively). All these considerations are readily applicable to the projection pruning of any product-lattice. 

Finally, we note that, as the result of the present complete pruning, we arrive at the set of admissible lattices and 
corresponding quasi-equilibria and equilibria. The question of which LB model can be supported by a particular 
sub-Maxwell lattice remains beyond the scope of this analysis. However, this is easily done upon studying the set 
of moments retained after the pruning. Note that, in general, the familiar single relaxation time lattice BGK model 
may be not sufficient, and more general kinetic models need to be addressed, such as the quasi-cquilibrium models 
[23l . [27l . |28| . [29| which make use of the quasi-equilibrium along with the equilibrium, or the multiple relaxation times 
(MRT) models (see, e. g., a paper by I. Ginzburg [3(| and references therein). As an illustration, a two-step quasi- 
cquilibrium model for incompressible flow is derived for the D3Q13 lattice in Appendix[El utilizing the above UniQuE 
quasi-equilibrium (|11|) . 

V. DISCUSSION 

Maxwell's derivation of the equilibrium distribution function in a gas predated Boltzmann's fundamental i7-theorem 
and the specification of the equilibrium as the minimum of H. Maxwell's argument was based on the independence of 
the equilibrium on the direction, resulting from the multiplication of the unidirectional equilibrium functions. Both 
approaches, Maxwell's and Boltzmann's, result in the same Gaussian equilibrium. 

In this paper, we followed closely the Maxwell's path, exploiting the symmetry of the product for the purpose 
of constructing LB models. The main result of the present theory is the constructive approach to better, Galilean 
invariant higher-order LB models. Here we summarize the construction of LB developed above, and make further 
comments on these findings. 

• Construction of any three-dimensional LB takes it origin in one dimension. For a given lattice, we consider the 
closure relation and derive the reference temperature via the matching condition. The reference temperature 
does not change in any further step of the construction. It is the characteristics of the one-dimensional lattice, 
and of all the lattices induced by the one-dimensional lattice in three dimensions (Maxwell and sub-Maxwell 
lattices). 

• Let us give another interpretation of the matching condition. The Maxwell moments arise from the Gaussian 
distribution ((T|). That means, they obey a recurrence relation which expresses the higher-order moments in 
terms of the two lower moments (the mean and the variance). This recurrence relation is well known and is 
not reproduced here. Important is that the moments of the Gaussian prolong: Once the first and the second 
moments are known, the rest of the moments are computed from the recurrence relation. Now, with a finite 
number of velocities Q (odd), we can reproduce first Q moments of the Gaussian (including normalization). 
However, this does not say anything yet whether or not the moment sequence will be prolonged. Because of 
the closure relation, such a prolongation is restricted to M(q) and M(q +1 ), the former is odd and was used 
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in the matching condition, the latter is even and leads to the same matching condition. Thus, the matching 
condition verifies a restricted Gaussian prolongation, it checks the moments which are not independent of the 
first Q moments (by the closure relation). But higher-order moments of the Gaussian are also dependent on 
the lower-order moments (through recurrence relation). So, the matching condition seeks consistence between 
the two different relations, the one is the closure relation (pertinent to the discreteness of the velocities), and 
the other pertinent to the Gaussian. This verification of the restricted prolongation is thus the verification 
of the restricted Gaussian feature for the given velocity set, and it reduces to the verification of the reference 
temperature, as it was done in sec. |TTJ In other words, important is not the reproducing of the Q moments of 
the Gaussian with Q populations (this can be done by any velocity set) but rather the prolongation property, 
which is the matching condition. 

• Transition to three dimensions begins with the construction of the Maxwell lattice, for which one defines UniQuE, 
the special quasi-equilibrium in the form of a product of unidirectional functions. UniQuE has remarkably simple 
moment relations (products of unidirectional moments), and reduces the analysis of the moment systems from 
Q 3 to 3Q — 2 dimensions. Construction of the equilibrium on the Maxwell lattice requires only the unidirectional 
information. 

• All other lattices are obtained as a pruning of the product lattice. The method of key-table reduces the 
problem of constructing the entropy function of the pruned lattice to analyzing a two-dimensional table and 
verifying consistency and solving linear problems. For large Q, this can be achieved with standard tools of 
linear programming (verification of consistency of linear systems). However, even the intuitive search for good 
sub-Maxwell lattices is possible with the key-table thanks to its relative simplicity. 

• Finally, the projection pruning is introduced as an extension of the key-table, which enables to derive UniQuE 
and Maxwcllization for pruned lattices. This requires only solving linear systems. Maxwellization on the pruned 
lattices is a promising approach to higher-order lattices due to a relative simplicity of construction. 

• Derivation of any lattice in any dimension begins with finding the reference temperature and the equilibrium 
at zero velocity (weights). After that, there are two options to continue, equilibration or Maxwcllization. 
The strong point about equilibration is that it is based on the entropy minimization, and stability theorems 
(Boltzmann's TJ-theorems) can be proved in that case for various LB realizations. However, in order to obtain 
the equilibrium on that route, one needs to solve a nonlinear minimization problem which, in most cases, can 
be only done within an approximation. On the other hand, in the Maxwellization approach, the corresponding 
equilibrium is constructed much easier, even for sub-Maxwell lattices it requires only solving linear systems. 

• We note that specific cases of UniQuE were used recently in order to construct quasi-equilibrium LB models 
with enhanced stability [ID, [3l[ , and to enhance Galilean invariance of LB models on standard lattices [l2| ■ 

• Finally, we point out that UniQuE represents an exact and systematic alternative to other closure procedures 
reported in literature, not necessarily in the LB context. For example, a moment-inversion algorithm was 
developed recently based on Cholcsky decomposition of the velocity covariance matrix and repeated applica- 
tion of one-dimensional quadrature for dilute gas-particle flows [32j. Even though it is well recognized that 
the moment-inversion problem admits exact solution in one dimension (e.g. by product-difference algorithm), 
defining the linear system used to solve for the weights and the abscissas in multi-dimensional case is still an 
open question. In particular, it is recognized [32Tj that the most suitable algorithms (Cholesky decomposition, 
method of eigenvectors etc) appears as problem dependent. In this context, UniQuE offers a simple and general 
framework to develop closure models starting from the analytical one-dimensional solution. Fixed abscissas used 
by Maxwell lattices do not represent a limit, since the closure relations in terms of the considered moments can 
be derived explicitly and implemented in functional form in the generalized hydrodynamic equations. Outcomes 
of this procedure are expected for granular flows, polydisperse liquid sprays undergoing droplet coalescence and 
evaporation and, more generally, aerosol dynamics [32|. These problems will be addressed in our future work. 
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APPENDIX A: HOW TO FIND CLOSURE RELATION AND VERIFY REFERENCE TEMPERATURE 

FOR A GIVEN VELOCITY SET 

For Q velocities (Q odd), one writes = aQ^v^ 2 + a (Q-4) l '(f) 4 + ' ' ' + a i w (i)i substitutes (Q— 1)/2 different non- 
zero values for the velocities and solves the linear system for the coefficients clq-2, ■ ■ ■ , ai- Once the latter are obtained, 
we use M ( M = b n T^ n ~ 1)/2 u + 0{u 3 ) {n odd) with b n = 1 x 3 x 5 • • • x n. Matching condition of linear in u terms results 

in the algebraic equation for the reference temperature, bgT^ — aQ^bq^T^ 3 ^ 2 — a^Q-^UQ^^T^ 5 ^ 2 — 
■ • • — ai =0. Positive roots (if they exist) define the reference temperature. If no positive roots are available, the 
corresponding lattice is ruled out of a further consideration. 



APPENDIX B: MAXWELL LATTICES AND ROOTS OF HERMITE POLYNOMIALS 

In Ref. [l[, it was argued that one-dimensional Maxwell lattices have ratios of the velocities that approximate the 
ratios of the roots of Hcrmite polynomials. We recover this argument here from the closure relation and the matching 
condition, considering the example of D1Q5. Without loss of generality, the one-dimensional velocities are set as 
V = {0, ±l,±r}, where r > 1. The closure relation then reads: v%-, = (1 + r 2 )v^ — r 2 v^y The matching condition 

results in the following quadratic equation for the reference temperature: 15T] 2 — 3(1 + r 2 )To + r 2 = 0. This equation 
has positive real- valued solutions if r > r* , where r* = \/t* with t* the larger root of another quadratic equation, 
3(1 + i) 2 - 20t = 0. From the latter we find t* = (7+2vT0)/3, and taking the root of it, we find r* = ( + Vf)/V3. 
This is nothing but the ratio between the two non-trivial roots of the 5-th order Hermite polynomial (these roots arc 

{0, ±a/5 - VlO, ±y/h + VlO}). Thus, we have recovered the argument of Ref. [H by a different consideration. 
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APPENDIX C: MAIN THEOREM ABOUT UNIQUE 

We here give the proof of the theorem of sec. IIIII which characterizes the UniQuE population ((4]) as the quasi- 
cquilibrium. We restore to expanded notation: For D = 3, the density is defined as 

M-)' (Cl) 

iev jev kev 

while the unidirectional moments are 

pM^ = E E E v l)hi,m > n=l,...,Q-l, 

i£V j£V keV 

p< ] = E E E »&> w n = i,...,Q-i, (C2) 

ievjev kev 

P M™ = E E E n = 1, . - 1, 

iey jev lev 

(n) 

The moment densities MA , a = x,y, z are termed unidirectional in order to reflect the fact that only the .T-component 
of the three-dimensional velocity vector vuj q = (v<i)i'V(j),'V(k)) participates in the definition of Mx , while only 

(n) 

the y-component participates in the definition of My , etc. Finally, we denote 

m x = {m^,...M q - 1) ), 

H^MM.-.Mf 1 )}, (C3) 

^ = {mw,..,mP»}. 

Theorem: Let the parameters M. a take their values in the positivity domain. A4 a S O, a = x, y, z. Then the 
minimizcr of the entropy function H , 

under the constraints (|C1[) and (|C2[) is given by the product-function ((4]). 
Proof: The extremum condition is written 



n— 1 n— 1 n— 1 



where A is the Lagrange multiplier corresponding to the density constraint (|C1[) . and are the Lagrange multipliers 



corresponding to the unidirectional moment constraints (|C2[) . This can be rewritten as 

f(i,j,k) = P X (*) Y U) Z (k), (C6) 

with 

X W = W {i) exp h- 1 ^ + E > 
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Substituting (|C6[) into the constraints (|C1|) and (|C2[) . the latter becomes 



(C8) 

- 1, (C9) 



(e*m) (E%) (e%)W 

(E^)*m) (E f 0) ] (E z w) = m ^ » = !>•■■, Q 

VieV / yev / Wv / 

E<) y w ] (E x w) (E z H = *4 n) . »=!,•■■, q 

v ieV / Viev / \kev / 

[E (E X (o) (E y W = M i n) ' n=l,...,Q-l. (Cll) 

Vfcev / Kiev ) \jev J 



-1, (CIO) 



Equation (|C8[) admits a solution (the normalization condition), 



which implies for the rest of the conditions, Eqs. I|C9[) . (|C10|) and (|C11I) . 

E X W = X > 

iev 

^^I (i) =MW, n = l,...,Q-l, 



E F o-) = 1 - 



E Z W = !' 

fcev 

5>f fe) %)=Mj">, »=1,...,Q-1. 



(C12) 



E X W = X > E F 0-) = !. E Z « = !' ( C13 ) 

iev jGV feev 



(C14) 



(C15) 



(C16) 



Now, each of the problems (|C14j) . (|C15[) and (|C16[) is equivalent to the one-dimensional problem solved in sec. ITU and 
which defines the one-dimensional functions (p^(A4), and thus the solution of each of these problems separately is 
given by the unidirectional quasi-equilibrium, viz. 



(C17) 



X(i) 


= <^(i)(Mi 1 ),. 








..,Mf-% 






.,Mp)) 



With (|C17[) and (|C13[) . we find a solution in the form The proof is completed by reminding that the minimum 
of a convex function under a set of linear constraints is unique. 
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APPENDIX D: D3Q27 AND D3Q19: UNIQUE, EQUILIBRATION AND MAXWELLIZATION 

Here we collect various populations for the D3Q27 Maxwell lattice and for the D3Q19 sub-Maxwell lattice mentioned 
in the paper. The list begins with the UniQuE fll} for the D3Q27: 



/(0,0,0) 


= p(l- 


u x )(i-u y ){i-n z ), 


/(*±1,0,0) 


1 fTT 

= 2(n, 


±U X )p(l - ILy)(l - I1 Z ), 


/(0,±1,0) 


1 

= ~ 2 P(1 


-IL x )(ILy±Uy){l -n z ), 


J(0,0,±1) 




-n x )(i-n y )(n 2 ± Uz ), 


/(±1,±1,0) 


= W 


E ± Ux )(n a ± % )(i-n z ), 


/(0,±1.±1) 


1 , 


-n x )(ii y ±u y )(ii z ±u z ), 


/(±1,0±1) 


= 


,±u x )(i-n y )(u z ±u z ), 


/(±1,±1±1) 


= w 


j ± U x )(jly i Uy)(Tl z i U Z ) 



(Dl) 



Note that, when setting II 2 = in the nine populations, /,* s, o)> /(*o ±i o)' ano - /(*±i ±i o) PIP - we obtain the 
UniQuE on the two-dimensional D2Q9 lattice: 

f* 00) =p(i-n x )(i-iiy), 
/(±i,o) = ^( n * ± «x)p(i - n y ), 

1 (D2) 

/(0,±1) = 2-PK 1 - nx)^ i Uy), 



f*±i,±i) = 4^( n z ± «.T)(n y ± u 



This two-dimensional UniQuE was used in [231 ] for a construction of a class of two relaxation times models with 
enhanced stability. 

Equilibration of (|D1[) is achieved upon substituting the equilibrium one-dimensional pressure. 

= \ (2^TT3^ - l) , (D3) 



3 



into dDl 



27 

pE 



/(o,o,o) = 7^P (2 - [2-<Jl + 31*2) ( 2 - + , 



/(±i,o,o) = ^fP (20 + 3u* - 1 ± 3 M:c J (2 - ^1 + 3u»J (2 - + 3 U 2 



2 



/fo.ibO) = (2 - (2^1 + 3 U 2 _ i ± 3% ) (2 - ^1 + 3^!) , 

2 

27 ; 



/fo.o.ii) = ^P (2 - v/1 + 3^1) (2 - v/l + (2^1 + 3^1 - 1 ± 3u, 



/(±i,±i,o) = (2^1 + 3u* - 1 ± 3u x ) (2^1 + 3^ - 1 ± 3 % ) (2 - Vl + 3^) , 



(D4) 



1 



/( E o,±i,±i) = (2 - v^+3^2) (2^1 + 3u 2 y -l±3u y ) [2^TT3^ z -l±3u z ), 
/(±i,o±i) = (2^1 + 3^2 - 1 ± 3u x ) (2 - ^1 + 3^) (2v/l + 3«2- 1 ± 3u 



/(±i,±i±i) = ^P ( 2 ^1 + 3«2 _ 1 ± 3u x ) (2y/l + 3ul-l±3uy) (2^1 + 3^ _ 1 ± 3 U; 
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Weights W s , corresponding to various shells s = 0,1,2,3 (sec (|9])), are numerical pre-factors in these expressions. 
Equilibrium (|D4p was derived in (33l | by a direct minimization of entropy in three dimensions. Positivity domain 
of (|D4|) (all populations arc non-ncgativc simultaneously) is a box with the edge 2 centered at the origin of the 
three-dimensional parameter space (u x , u y ,u z ): ^d 3 q 2 7 = {u : \u a \ < l,a = x,y, z}. 

Maxwellization of (|D1|) is found upon a substitution into (|D1|) the Maxwell expression for diagonal components of 
the pressure tensor at unit density, 



C = ±(l + 3^), (D5) 



which gives explicitly 



21' V 2 x / \ 2 
2 

1,0) ~ 27/ 



/(±i,o,o) = ^P (1 ± 3u a + 3u£) ( 1 _ luA (1 -\u 3 B ), 



/(o, ±ll o) = J^P (l - f «*) (1 ± 3% + H) i 1 - \< 

— * „ ( 1 3 „,2 \ 



/(o,o,±i) = ^P ( 1 - 5 «2 ) ( 1 - 5 «S ) C 1 ± 3 ^ + 3u 



"z J 1 



/(±i,±i,o) = ^P (1 ± + 3t£) (1 ± 3 % + 3u*) f 1 - -u: 



_ J. „ /1 _i_ q„, 1 q„,2\ /1 _i_ Q „, , q„,2\ ( ! 3 „,; 
5- 

/(o,±i.±D = ^P ( 1 - ) C 1 ± 3 «y + 3u ') (! ± 3w * + 3u ') . 



(D6) 



/(±i,o ± D = ^p (1 ± 3 «* + 3 «') f 1 - §«2 ) (1 ± 3 ^ + 3 " 2 



/(±i,±i±i) = 2leP (1 ± 3 «x + 3u») (1 ± 3u y + 3u 2 y ) (1 ± 3u z + 3u 2 z ) . 

Positivity domain of (|D6|) is the box with the edge 2y2/3: ^D3Q27 = i u ■ \ u a\ < \/2/3, a ~ x iV^ z }- 
For the D3Q19 sub-Maxwell lattice, the analog of UniQuE constructed by projection pruning (jTTJ) is: 

/(o,o,o) = p(! - n - - n y - n z + + n 9 n 2 + n x n z ), 
/(±i,o,o) = ^p(! - n y - n ^)( n x ± 

/(0,±1,0) = \p( 1 - Ha; - n 2 )(n y ± My), 

/(o,o,±i) = ^p(! - n x - n y )(n 2 ± ( D7 ) 
/*±i,±i,o) = \pi^x ± ?^)(n y ± Uy), 

/(*o,±i,±i) = ^P( n y ± u y)( n z ± «*), 

/(*±i,o±i) = ^p( n * ± «x)(n 2 ± u z ). 

It is easy to verify by a direct computation that the moments of the populations (|D7j) satisfy the relations given by 
the second row of Tabic ITTT1 Note that, when setting H z = in the nine populations, f? j, f* ±1 j, /* ±1 ), and 
f(±l ±1 0) 4EZl) we again obtain the UniQuE on the two-dimensional D2Q9 lattice (|D2[) . Maxwellization (Ti"3"|) of (|E>7|) 



18 



is achieved upon substitution of (|D5 



/(o,o,o) = \p[ 1 - ( u l +u l + U D + 3 ( u l u l + u l u l + U K)] > 
/(±i,o,o) = ^P [1 - Kul + ul)] (1 ± 3u a + 3nl 



"X ) 1 



/(o,±i,o) = t 1 - 3(t£ + u»)] (1 ± toy + toft > 

/(o.o I± i) = [1 - 3 (<4 + C 1 ± 3 «- + 3 «') > ( D8 ) 

/(±i,±i,o) = C 1 ± 3m * + 3u ?) (l ± 3w y + 3u*) , 

/(o 1± i,±i) = (1 ± 3 % + 3 «y) (1 ± 3u z + iu 2 z ) , 



/(±i,o±i) = ™P (1 ± 3m * + 3 <4) (1 ± 3m, + 3u^ 



1 

36 r 

Positivity domain of (|D8[) is the intersection of three cylinders, C x = {u : u y + uj. < |}, C y = {u : u x + u 2 z < i} and 

= {u : u 2 + u y < J,}: ^D3Qi9 = C x f]Cy f] ^z- Since ^D3Qi9 i s included in a box with the edge 2/v / 3, we have 
the following inclusion relations between the positivity domains: 

^D3Q19 C ^D3Q27 C ^D3Q27- (D9) 

Although the positivity domain shrinks when proceeding from the Maxwell to the sub-Maxwell lattice, all the three 
equilibria are well consistent with the low Mach number restriction to these models, |« a | <C 1/V3- Functions (|D7|) 
and (ID8I) are used in [3l| for the construction of a three-dimensional two relaxation time LB model. 



APPENDIX E: QUASI-EQUILIBRIUM D3Q13 MODEL 

The D3Q13 is the sub-Maxwell lattice of the D3Q27 with the smallest number of velocities capable of retaining 
the pressure tensor. The peculiarity of the D3Q13 as compared to the other lattices (the Maxwell D3Q27 and the 
sub-Maxwell D3Q15 and D3Q19 lattices) is in the third-order moment tensor Q a B-y- Indeed, the D3Q27, D3Q15 and 
D3Q19 lattices all recover the isotropic linear part of the equilibrium function Q„g^ in the form 

<2a/3 7 = 7^(u«<Vy + u o5 a ~f + Ury5 a p) + 0(u 3 ), (El) 

which corresponds to the linear in u piece of the correct Maxwell moment relation at the reference temperature 
To = 1/3. Terms of order 0(u 3 ) are different for each of the D3Q27, D3Q15 or D3Q19 lattices but their effect is 
negligible at low Mach numbers. On the contrary, the corresponding expression for D3Q13 is not isotropic even at 
the linear order: 

u a , 

^u a + 0(u 3 ), (E2) 
0(u 3 ). 

Note that the factor 1/2 instead of 1/3 in the off-diagonal terms Q^pp (|E2[) is inconsistent with the correct Maxwell 
relation (|E1|) (in other words, the diagonal terms Q^ aa in (|E2j) correspond to the correct reference temperature 
To = 1/3 whereas the off-diagonal terms Q^bb correspond to a different "temperature" 1/2). Thus, the D3Q13 lattice 
is less isotropic than any of the other sub-Maxwell lattices (D3Q15 or D3Q19) of the Maxwell D3Q27 lattice. This 
peculiarity precludes developing the standard LBGK model on the D3Q13 lattice, as was first noticed in [2(| upon a 
different consideration. 

Utilizing the concept of UniQuE, we shall now derive a simple BGK-like model with two relaxation times which 
recovers the incompressible Navier-Stokes equations on the D3Q13 lattice. For that, we use a generic pattern of 
quasi-equilibrium kinetic equations with a two-step relaxation mechanism (23l . [27l . [28l . [23 . l3~fl ] , 

d t f + v-Vf = --(f-n--(r-f M ), (E3) 



Qa 



M 

aa 



^aBB 
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where the first term in the right hand side describes a relaxation to the UniQuE state /* (with a rate ti), and 
the second term represents a relaxation from the UniQuE to the equilibrium (with a rate r 2 ). A rationale behind 
using a two-step quasi-equilibrium model (|E3j) in the present context is the following: The two steps of relaxation 
"adjust" separately the off-diagonal and the diagonal components of the nonequilibrium pressure tensor (see below) 
and will be tailored in such a way as to recover isotropy in the low Mach number limit (see also Refs. j23l [29[ for 
the application of this type of models in various other context). For the present case, we choose the UniQuE of the 
projection pruning (see Tab. IIII[) , and the equilibrium as the Maxwcllization thereof. Specifically, the UniQuE (JTTJ) is 
found by an inversion of the moment relations given in the D3Q13 row of Tab. IIIII It proves convenient to restore a 
notation Tl aa = Tl a for the diagonal components of the pressure tensor at unit density: 



f (0,0,0) 


= p fl - l;(ttxx + n ra 


+ n 2Z ) 


)■ 








f(a,\,0) 


= -p[(n xx + au x )(l + 


XUy) + 




\-\Uy){l- 


1- <TU X ) - U zz (l + <TU X - 


f \Uy)] , 


f(a,0,X) 


= -p[(Il xx + au x )(l + 


Xu z ) + 


(n„ - 


h Xu z )(l - 


- au x ) - Tlyy ( 1 + au x - 


h Xu z )} , 


f(0,v,X) 


= zPii^yy + < 7U y)( 1 + 


Xu z ) + 


(n 22 H 


-Au 2 )(H 


- auy) - E xx {\ + au y - 


- Xu z )} , 



(E4) 



where a £ { — 1, 1} and S £ { — 1, 1}, and the Maxwcllization (|13|) is achieved upon substituting 11^ = 1/3 + u 2 a (|D5j) 
into the above expression (|E4[) . 

It can be shown that, if the relaxation times T\ and T2 are chosen as 

n = r, r 2 = 2r, t > 0, (E5) 

then, under the diffusive sca ling at low Mach number (dt — ► e 2 dt, d x — > ed x , u = eu^\ p = 1 + e 2 p < - 2 \ where e is 
the Mach number, see, e. g. [34. l35l [36l. [37l. [38| ) . kinetic equation (|E3|) reduces to the incompressible Navier-Stokes 
equation, 

V-m (1) =0, (E6) 
+ • Vtt (1) + Vp - j/A« (1) = 0, (E7) 

where A = d 2 + d 2 + d 2 is Laplace operator, p = p 1 - 2 ^ /3 is the hydrodynamic pressure defined by the solenoidal 
(incompressibility) condition (|E6j) . and v is the kinematic viscosity given by the formula, 

v = t/2. (E8) 



Note that the kinetic model (|E3[) is realizable under the condition (|E5|) : relaxation towards the quasi-equilibrium is 
faster than the relaxation from the quasi-equilibrium to the equilibrium (ti < t 2 ). 

The simplest way to prove this statement is to consider a closed moment system equivalent to the kinetic model 
(IE3D for the moments 



pM = {p, pu x , pu yi pu z , pIL xx , pUyy, pTL zz , pH xy , pRyz, P^xz, pT x , pT y , pT z }, (E9) 
where the three independent third-order moments T a arc defined as 

P T * = ( v (i)( v U)- v2 (k))f( t ,3,k)), 

P T v = \ v (j)( v h - v !k))f^,3.k)) , (E10) 
pT z = (u(fc)(W(i) ~ rfj))f(ij,k)) • 
The moment system equivalent to the kinetic equation (|E3[) reads 

e 2 d t p + ed x (pu x ) + ed y (pu y ) + ed z (pu z ) = 0, (Ell) 



e 2 d t (pu x ) + ed x (pU xx ) + ed y (pll x y) + ed z (pll xz ) = 0, 

e 2 dt(pu y ) + ed x (pU xy ) + edy(pUy y ) + ed^Ey,) = 0, (E12) 
e 2 d t {pu z ) + ed x { P n xz ) + edy(pUy Z ) + ed z (pU zz ) = 0, 
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e d t (pUxx) + edx(pu x ) + -ed y (p(u y + T y )) + -ed z {p(u z + T z )) = ~~P [ ^xx - ( ^ + ' 
e 2 d t {pU yy ) + ed y (pu y ) + ^ed x (p(u x + T x )) + \d z {p{u y + T y )) = ~p ( 
e 2 d t (pU zz ) + ed z (pu z ) + \ed y (p{u y + T y )) + \ed z {p(u z + T z )) = ~p (n zz - (l + u \ 

e 2 dt{pH X y) + \idx{p{u y + T y )) + l;edy(p(u x + T x )) = -— p(U xy - u x u y ), 

Z Z T\ 

e 2 d t (pU yz ) + ^ed y {p(u z + T z )) + \ed z {p{u y + T y )) = ~p{Tl yz - u y u z ), (E14) 
e 2 d t ( P n xz ) + \ed x {p(u z + T z )) + -ed z (p(u x + T x )) = -—p{U xz - u x u z ), 

z z n 

e 2 dt(pT x ) + ed x (p{n yy - H„)) + ed y {pll yz ) - ed z (pIL xz ) = - -p(T x - u x {Il yy - U zz )) 

n 

- —pu x (Ji yy - n zz -ul + u 2 z ), 

e 2 d t {pT y ) + edy(p(U xx - iy) + ed x (pU X y) - ed z {pH yz ) = —p(T y - u y (IL xx - H zz )) 

n 



(E15) 

pUy(jlxx - H ZZ ~U 2 + U z ), 

e 2 d t (pT z ) + ed z (p{n xx - U yy )) + ed x (pll xz ) - ed y {pU yz ) = - —p{T z - u z (H xx - U yy )) 

- — pu z {\\ X x - n TO -ul + ul), 

where we have explicitly introduced the diffusion scaling. Substituting p = l + e 2 p( 2 > and u = eu^ into the continuity 
equation (|Ellj) . we find at the first non-trivial order: 



= 0, (E16) 



where summation convention is applied. Next, from the relaxation equations for the components of the pressure 
tensor, Eqs. (|E13)) and (|El4]) . it follows that 

pU aP = l -5 aP + e 2 Q<W 2) + + * 2 <? {2) + °( £3 )> (E17) 

where the term n^^ 2 -* is the non-equilibrium (viscous) part of the pressure tensor which is not yet defined. Substi- 
tuting the latter expression into the momentum equation (IE12I) yields 



(d t u^ + d a (pW/3) + ufdpuW + ^n^ q(2) ) = 0, (E18) 



where we have made use of the solenoidal condition (|E16[) . What remains is to derive the nonequilibrium part fl^ q . 
For that, let us consider again the moment equations for the components of the pressure tensor. These give: 



e 2 \(dxUyV+d y u^)=-±e 2 n^ 2 \ (E19) 



for the off-diagonal component IL xy and similarly for the rest of the off-diagonal components, Eq. (|E14|1 , and 

W> + \ (dy^ ] + = -y/K^ 2 \ (E20) 



r 2 



for the diagonal component H xx and similarly to other diagonal components, Eq. (|E13|) . Note that, when deriving 
the above results, we have used the property of the third-order moments, T a — 0(e 3 ), which follows from the right 
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hand side of the moment equations for T a (Eq. (|E15p ). Using once again the solenoidal condition (|E16|) . Eq. (|E20|) 
can be rewritten: 

\ e *(d x uV + d x uV)=-^Ul^\ (E21) 
Thus, by choosing the relaxation times as T\ = r, r 2 = 2r, the noncquilibrium pressure tensor becomes isotropic: 

T neq(2) 



rr 



~(d a uf +d p u^). (E22) 



Substituting (|E22[) into the momentum equation (|E18[) concludes the derivation of the incompressible Navier-Stokcs 
equations (|E7|) from the quasi-equilibrium kinetic model (|E3[) . 

Finally, it is strai ghtf orward to derive a lattice Boltzmann scheme for the kinetic equation (|E3j) following a general 
method of Refs. 23|, 2fJ|3lj]: Kinetic equation (|E3|) is integrated in time from t to t + St along characteristics, and the 
time integral of the right hand side, J = — — (/ — /*) — — / )> is evaluated by trapezoidal rule to get 

f(x + vSt, t + St)- f(x, t) = ^ J(f(x + vSt, t + St)) + j J{f{x, t)). (E23) 
In order to avoid implicit computations in the latter expression, let us apply the following variable transform [Til. l29j|: 

f^9 = f-jJ(f), (E24) 
to Eq. (|E23|) . which yields after taking into account (|E5|) : 

g(x + vSt, t + St) = (l- uj)g(x, t) + | [f M (p, u) + f* (p, «, K a )] , (E25) 



where 

2St 



2t + St ' 

9 = P(9), (E26) 
u = u(g), 

Ka = [^rU aa (g) + 5tn™ a (g)] . 

The scheme (|E25|) becomes the LB scheme if the time step St is matched with the lattice. A different MKT LB 
equation for the D3Q13 lattice was suggested in [2(| . 



